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Abstract 


A  method  for  modeling  images  of  natural  terrain  is  developed  and  applied 
to  the  segmentation  of  aerial  photographic  data.  An  underlying  stochastic 
structure  based  on  linear  filtering  concepts  provides  a  means  of  modeling 
the  terrain  in  local  areas  of  the  image.  Superimposed  on  this  is  a  Markov 
random  field  that  describes  transitions  from  regions  of  one  terrain  type 
to  another.  Maximum  likelihood  and  maximum  a  posteriori  estimation  is 
applied  to  estimate  regions  of  similar  terrain.  Results  of  application  to 
digitized  aerial  photographs  of  a  rural  area  are  presented  and  discussed. 
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I-  INTRODUCTION 


The  segmentation  of  images  is  a  vital  part  of  image  analysis  for  a  host 
of  applications.  Where  specific  models  for  images  are  lacking,  one  is 
forced  to  base  the  analysis  on  heuristic  arguments  and  heuristically 
motivated  features  and  performance  can  be  evaluated  only  on  an  empirical 
basis.  Much  of  the  work  in  image  analysis  has  fallen  into  this  category. 

When  one  is  able  to  restrict  the  class  of  images  so  that  specific  models  can 
be  developed,  the  resulting  analysis  algorithms  are  placed  on  a  much  firmer 
foundation  and  their  performance  can  be  evaluated  from  both  analytical  and 
empirical  considerations.  We  have  taken  the  latter  approach  in  this  report. 

Many  images  of  natural  terrain  can  be  effectively  characterized  as  two- 
dimensional  (2-D)  random  processes.  The  segmentation  of  these  images  into 
areas  of  known  types  is  important  for  a  variety  of  applications  including 
military  surveillance  and  ’■econnaissance ,  crop  and  land  use  data  collection, 
and  cartography.  In  addition,  many  images  other  than  those  of  terrain 
contain  areas  of  texture  that  can  be  well  modeled  by  random  processes. 

This  report  describes  a  model  for  images  based  on  stochastic  linear  filtering 
concepts.  Using  this  model,  one  can  develop  the  probability  density 
functions  for  the  image  data.  Segmentation  of  the  image  is  then  treated  as 
an  estimation  problem  and  the  resulting  algorithms  have  a  clear  intuitive 
interpretation. 

Earlier  results  of  this  research  as  applied  to  texture  data  were  reported 
in  Reference  1.  This  report  extends  those  results  and  demonstrates  their 
application  to  terrain  image  data.  While  no  claim  of  universality  is 
made,  the  algorithms  have  been  effective  in  segmenting  aerial  photographs 
of  rural  terrain  and  provide  encouragement  for  continued  research  along  these 
lines. 

In  the  following  section,  we  postulate  a  general  class  of  linear  filtering 
models  and  then  specialize  the  results  somewhat  for  the  work  covered  in  this 
report.  Next,  we  consider  segmentation  as  a  region  estimation  problem  and 
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develop  maximum  likelihood  (ML)  and  maximum  a  posteriori  (MAP)  estimates  for 
the  regions  (segments).  Finally,  vt  present  some  results  of  application  of 
the  ML  and  MAP  algorithms  to  aerial  photographic  data  and  discuss  and  compare 
the  two  approaches. 
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II.  LINEAR  FILTERING  MODELS  FOR  IMAGES 
A.  In  this  section,  we  consider  models  for  an  image  of  the  form 

F'(n,m)  =  a  F'(n-i,m-j)  +  ]C  b  W(n-i,m-j)  (la) 

a  3  e  3 

(i,j)¥(0,0) 

F (n,m)  =  F’(n,m)  +  G  (lb) 

ne[0,N-l] 
me [0,M-1] 

where  a  and  B  are  finite-extent  masks  covering  the  filtered  points,  W(n,m) 
are  a  set  of  independent  identically  distributed  zero  mean  random  variables, 
and  G  is  a  constant  representing  the  mean  value  of  the  image.  In  general, 

G  could  be  made  a  function  of  the  image  coordinates  (n,m)  but  this  implies  a 
form  of  non-stationarity  and  an  associated  registration  problem.  Since  the 
images  to  be  dealt  with  are  (at  least  locally)  stationary,  we  shall  not  make 
that  generalization. 

In  the  following,  it  will  be  assumed  that  a  and  S  are  chosen  so  that  F' 
is  recursively  computable  [2,3]  given  W  and  vice  versa.  This  will  lead  to 
simplifications  in  the  analysis  to  follow. 

If  the  vectors  f_,  w,  and  £  represent  an  ordered  set  of  the  corresponding 
array  points  (derived  for  example  by  scanning  rows  of  the  image),  then 
Eq.  (1)  can  be  written  in  a  matrix  formulation  as 

A(f-g)  =  B  w  +  x  (2) 

where  A  and  B  are  matrices  whose  non-zero  elements  are  derived  from  the  terms 
a^j  and  b^  in  Eq.  (1)  and  represents  a  set  of  boundary  conditions  imposed 
by  the  finite  extent  of  the  image.  In  general,  these  boundary  conditions 
should  be  modelled  as  random  variables  or  the  parameters  of  the  model  should 
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be  made  space-varying  over  the  region  where  the  spatial  masks  a  and  6  are 
not  fully  contained  within  the  boundaries  of  the  image  (see  Ref.  4  for  a 
discussion  of  this  problem  with  respect  to  1-D  signals).  The  inclusion  of 
these  effects,  however,  is  of  little  significance  in  most  problems  of 
practical  interest  and  merely  serves  to  complicate  the  analysis.  Thus,  it 
is  convenient  to  set  the  boundary  conditions  _y  to  _0  and  accept  the  approxi¬ 
mation  that  results.  Since  the  terms  W(n,m)  are  independent,  one  can 
solve  Eq.  (2)  for  w  and  express  the  multivariate  probability  density  function 
for  the  image  as 

Pf<£)  =  . ~i~'7  P_.(B  1a(£-^)) 

-  |A  B|  - 


(N-l.M-l) 

n  Pw(E(n,m)) 

(n,ra)=(0,0) 


(3) 


where  the  notation  E(n,m)  is  used  to  represent  the  ordered  components  of 
the  vector  B  ^A(f-j>)  and  where  we  have  used  the  fact  that  since  both  A  and 
B  correspond  to  recursively  computable  masks,  their  determinants  are  equal 
to  one.  Further,  since  mask  6  has  a  recursive  form,  E(n,m)  can  be  computed 
from 


E(n,m)  =  £  b  E(n-i,m-j )  a  F’(n-i,m-j)  (4) 

B  J  ct  1J 

(i,J)*(0,0) 
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Equation  (3)  provides  a  formula  for  the  probability  density  function  of 

image  F  in  terms  of  the  probability  density  function  for  the  driving 

variables  W(n,m)  and  is  key  to  the  estimation  theoretic  approaches  that  are 

taken  in  this  report.  Later,  explicit  forms  for  the  filter  and  the  density 

p  will  be  assumed  and  more  specific  results  based  on  these  forms  will  be 
w 

derived.  However,  the  methods  just  discussed  are  capable  of  representing 
only  images  with  homogeneous  properties  such  as  an  image  of  a  single  type 
of  terrain*.  The  images  to  be  considered  for  segmentation  contain  a  number 
of  different  homogeneous  regions.  Thus,  it  is  necessary  to  add  another 
level  of  structure  to  the  model  to  represent  these  more  typical  images.  This 
topic  will  be  discussed  next. 

B.  Consider  an  image  such  as  that  depicted  in  Figure  1  consisting  of 
regions  R^.R^, . . . ,Rp.  It  will  be  assumed  that  wl thin  each  region,  the 
image  can  be  represented  by  a  model  of  the  type  discussed  in  Part  A.  Denote 
the  probability  density  function  for  the  driving  terms  in  a  region  R^  of  type 
k.  by  p  (.).  Note  that  there  may  be  more  than  one  region  of  a  given  type. 

1  i 

Then,  ignoring  boundary  effects,  we  can  express  the  probability  density  for 
the  image,  given  the  regions  R^ , . . . ,R^  as 


p  (F  |  R  i,...,Rp)  =/n  p  (E  (n,m)  )\  •(  II  p  (E  (n,m))\*. 

\(n,m)eR1  kl  kl  /  \n,m)CR2  k2  k2  / 


/n  Pk  (n,m))j 

\n,m)cRp  >  p  / 


(5) 


*  In  this  report,  we  shall  refer  to  images  whose  statistical  characteristics 
are  spatially  invariant  as  "stationary"  and  use  the  term  "homogeneous"  only 
in  a  loose  intuitive  sense.  Thus,  the  word  homogeneous  should  not  be  con¬ 
fused  with  the  mathematical  property  used  in  the  literature  on  random  field 


where  the  subscript  k^  on  the  residual  terms  E(n,m)  indicates  that  these  terms 
wt  ,-e  computed  using  the  filter  of  type  k^.  Equation  (5)  or  its  logarithm 
can  be  used  as  a  likelihood  function  to  generate  maximum  likelihood  estimates 
for  the  regions.  This  procedure  will  be  detailed  in  section  III.  However, 
maximum  likelihood  estimation  implicitly  assumes  that  all  region  configurations 
are  equally  likely  and  this  is  seldom  the  case  for  most  images.  In  reality, 
regions  for  a  given  class  of  images  may  be  more  1 ikelv  to  be  of  some  approxi¬ 
mate  size  and  density  and  it  is  not  unusual  to  have  some  rough  prior  infor¬ 
mation  about  the  occurrence  of  regions  within  the  image.  When  such  prior 
information  is  available,  a  better  procedure  is  to  use  some  form  of  Bayes 
estimation.  This  can  be  approached  as  follows. 

Let  the  occurrence  of  regions  within  an  image  be  modeled  through  region 
transition  statistics.  In  particular,  define  the  "state"  s(n,m)  of  a  point 
(n,m)  as  the  region  type  to  which  that  point  has  been  assigned.  Thus,  if 
there  are  K  region  types,  the  state  may  take  on  values  in  the  set  of  integers 
il,2,...,K}.  Next,  assume  that  the  state  of  a  point  is  stochastically 
dependent  on  some  adjacent  set  of  states  in  a  support  region  ^  (see  Fig.  2). 
In  particular,  following  Kaufman  et  al» [5] ,  we  shall  assume  for  the  moment 
that  the  states  form  a  Markov  chain  specified  by  transition  probabilities 
Pr[s(n,m)|Sn  mJ*.  The  Markov  chain  imposes  constraints  on  the  shape  of  the 
state  support  region.  For  example,  the  symmetric  region  indicated  by  the 
dashed  lines  in  Figure  2  is  not  allowed.  We  shall  relax  the  constraint  on 
Sn,m  later  bF  allowing  the  probabilistic  structure  of  the  states  to  have 
other  Markov  dependencies.  However,  the  Markov  chain  assumption  will  permit 
development  of  the  estimation  equations  for  the  state  assignments  in  a 
straightforward  manner.  Since  the  set  of  all  possible  state  assignments  is 
one  to  one  with  the  set  of  all  possible  divisions  of  the  image  into  regions, 
we  shall  have  reached  our  objective. 

*  Our  notation  is  convenient  but  may  be  somewhat  confusing.  A  proper  (but 

more  cumbersome)  notation  would  be  Pr  [s  (n  ,itO=SqI  m=Sgl  where  s^  is  an 

integer  in  the  set  {1,2, ...,K}  and  SQ  is  composed  of  integers  in  the  same 

set.  Thus,  when  we  write  P  r  [  1 1 S  ]  as  in  F.q  .  (14),  the  reader  should 

n,m 

interpret  this  as  Pr [s (n ,m)=l | m=SQ^' 
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Let  S  represent  a  particular  set  of  state  assignments  for  all  points  in 

an  image.  Our  goal  is  to  evaluate  Pr[S],  If  the  image  is  N  by  N  pixels  in 

2 

size,  then  let  the  pixels  be  numbered  from  1  to  N  according  to  the  ordering 
imposed  by  the  Harkov  chain.  Then  we  can  write 

Pr[S  ]  =  Pr[sN2  ,  SN2_1> . . .  .sj 

=  pr[sN2lsN2_i,...,Si]*PrtsN2_1lsN2_2,...,sil--*  Pr[s2|s1].Prls1] 


(6) 

Let  be  the  support  region  for  s^.  Then  when  is  completely  within  the 
boundary  of  the  image,  we  have  Pr  [s^  |  s^_^  , . . .  ,s^  ]  =  PrO^IS^].  When  is 
not  completely  within  the  image  boundary,  i.e.,  when  s^  does  not  have  full 
support,  we  shall  approximate  Pr [s^ [ s^_^ , . . . .s^ ]  by  a  term  of  the  same  form 
where  zeros  are  inserted  for  the  states  not  present.  Thus,  Eq.  (6)  can  be 
expressed  as 


N‘ 

Pr  [S  ]  =  2  Pr[s  |S,  ]  =  II  Pr[s(n,m)  |  Sn  J  (7) 

k=l  k  k  (n,m)  n,m 

in  image 


Equations  (5)  and  (7)  can  be  combined  to  develop  conditions  for  the  maximum 
a  posteriori  estimates  of  the  state  assignments  which  define  the  regions. 
This  is  discussed  in  the  next  section. 
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III.  SEGMENTATION  AS  A  REGION  ESTIMATION  PROBLEM 

The  basic  concepts  for  approaching  segmentation  as  a  region  estimation 
problem  were  developed  in  the  previous  section.  In  this  section,  we  will 
restrict  attention  to  a  particular  class  of  linear  filtering  models  and 
develop  some  more  specific  results. 

The  assumptions  that  are  made  in  the  following  are  that  the  model  of 
Eq.  (1)  is  purely  autoregressive  so  that  b  =  6^  (the  Kronecker  delta), 
that  the  mask  a  has  a  quarter  plane  causal  form,  i.e.,  a  is  non-zero  only 
for  ie[0,I]  and  je[0,J]  and  that  the  W(n,m)  are  Gaussian  with  density  function 


1  ,  w  N 

rr  exp 

/2ira  a 


For  the  causal  autoregressive  model,  the  filter  coefficients  can  be 
efficiently  computed  using  the  algorithm  described  in  Reference  6. 

Taking  minus  twice  the  log  of  Eq.  (5)  and  applying  Eq .  (8)  we  obtain 


V'  r  Rk  ^n,n^  2  "I 

-2  in  p(F|R1,R2>  . . .  ,Rp)  =  1 _  +  incs^  +  ... 

R1 L  o, 


E  (n,m) 

+  L  I  „  2 

R_  I  2  +  f  nCTk 


1  2 

-  N  <!n2Tr 


i,  E^  (n,m) 

-  >  )  K-i  2  2 

,  ,  ,  v  „  — 1 -  +  in  o,  -N  In  2 it  (9) 

i=l  (n,m)e^  Q  2  ^ 

L  ki 
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For  the  maximum  likelihood  procedure,  the  number  of  regions  P  and  the  regions 
themselves  are  assumed  to  be  deterministic  parameters  of  the  density  function. 
An  ML  estimate  for  these  parameters  is  obtained  by  choosing  values  that  maxi¬ 
mize  Eq.  (5)  or,  equivalently,  minimize  Eq.  (9).  A  moment's  reflection 
indicates  that  Eq.  (9)  is  minimized  if  every  point  (n,m)  in  the  image  is 
assigned  to  a  region  R^  or  type  k^  such  that  the  term  in  brackets  is  minimum. 
For  the  case  of  two  region  types  (but  possibly  multiple  regions),  we  are  led 
to  a  segmentation  rule  of  the  form 


E^n.m) 


Q 


> 

1  © 


E2(n,m) 


+  2na. 


(10) 


where  the  number  above  or  below  the  inequality  indicates  the  region  type  to 
which  the  point  (n,m)  will  be  assigned  if  the  inequality  holds. 

Since  the  ML  method  leads  to  a  decision  rule  that  assigns  points 
to  region  types  without  regard  to  the  assignment  of  adjacent  points,  one 
might  expect  this  segmentation  rule  to  produce  a  number  of  false  assignments 
leading  to  a  somewhat  "spotty"  result.  The  examples  in  the  next  section  show 
that  this  is  in  fact  the  case.  The  main  advantage  of  the  ML  method  is  that 
it  is  quick  and  easy  to  apply.  The  spottiness  of  the  result  can  some¬ 
times  be  removed  by  lowpass  filtering. 

An  improved  form  of  region  estimation  is  Bayes  estimation  where  the 
regions  are  treated  not  as  unknown  parameters,  but  as  random  quantities 
whose  statistical  properties  are  described  by  the  Markov  transition  model 
developed  in  the  previous  section.  In  particular,  we  will  consider  maximum 
£  posteriori  (MAP)  estimation  where  we  maximize  the  probability  of  a  given 
set  of  regions  conditioned  on  our  observation  of  the  image.  From  Bayes  rule, 
the  a  posteriori  probability  can  be  written  as 


Pr[R1,R2,...,Rp|F] 


p(F|R1,R2,...Rp)Pr[R1,R2 

pTf) 


(ID 
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Since  the  denominator  is  not  a  function  of  the  regions,  it  is  only  necessary 
to  be  concerned  with  the  terms  in  the  numeratin'.  Observe,  first,  that  since 
the  set  of  all  possible  region  assignments  for  an  image  is  one-to-one  with 
the  set  of  all  possible  state  assignments  for  the  pixels,  the  second  term 
in  the  numerator  of  Eq.  (11)  is  given  by  Eq .  (7).  Next,  observe  that  the 
double  summation  in  Eq.  (9)  can  be  written  as 


all  (n.oi)  02  S 

s(n,m) 


Now,  taking  minus  twice  the  log  of  Eq.  (11),  using  Eqs.  (7),  (9)  and 
(12),  and  eliminating  the  constant  terms,  we  find  the  MAP  estimate  requires 
that  the  pixel  states  be  chosen  to  minimize 


£ 

all  (n,m) 


E  (  \ (n,m)  „ 

’  +  in a  ,  v  -2  in  Pr[s(n,m)|s  )1 

- 2 -  s  (n,m)  *  1  n,m  ‘ 

°s(n,m) 


This  requires  that  each  term  in  the  sum  be  separately  minimized.  In  the 
case  of  two  region  types,  this  reduces  to  the  conditions 


E,(n,m)  vU 

- = - +  inoj  -2«nPr(l|Sn  J  < 

2  1  n,m  > 

°1 


E2(n,m) 


+  ino  -2in  Pr T 2  j S  ] 
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Since  the  states  are  assumed  to  fcrm  a  Markov  chain,  Eqs.  (13)  or  (14)  can 
be  solved  for  the  states  recursively  according  to  the  ordering  of  the  chain. 

The  segmentation  algorithms  as  thus  far  developed  have  the  characteristic 
that  all  of  the  processing  of  the  image  to  effect  the  segmentation  is  direc¬ 
tional.  Since  many  terrain  images  have  no  preferred  direction,  it  is  desir¬ 
able  to  perform  operations  on  the  image  which  are  omni-directional.  Unfor¬ 
tunately,  2-D  filters  whose  outputs  are  recursively  computable  do  not  allow 
for  this  symmetry.  However,  under  certain  conditions,  it  is  possible  to 
allow  the  Markov  transition  probabilities  to  be  defined  for  a  symmetric 
state  support  region.  Such  models  have  been  studied  in  the  context  of  Markov 
Random  Fields  [7-9], 

In  the  case  of  a  symmetric  state  support  region,  F.qs.  (13)  or  (14) 

can  no  longer  be  solved  recursively.  Solution  bv  direct  methods  for  any 

2 

reasonably  sized  images  appears  to  be  hopeless  since  all  of  the  N  states 
are  coupled  through  nonlinear  equations.  However,  the  equations  can  be 
solved  by  iteration  as  follows.  An  initial  set  of  states,  usually  the  maxi¬ 
mum  likelihood  states,  are  assigned  to  the  image  points.  The  Markov  trans¬ 
ition  probabilities  in  Eqs.  (13)  or  (14)  are  then  based  on  these  state 
assignments  and  held  fixed  while  the  equations  are  evaluated  to  determine 
a  new  set  of  states.  These  new  state  estimates  are  now  inserted  into  the 
equations  and  the  procedure  is  repeated  in  an  iterative  manner.  If  a  stage 
of  iteration  is  reached  where  the  state  assignments  no  longer  change,  then 
the  MAP  estimate  has  been  found.  If  such  a  stage  is  not  reached,  this  mav 
be  an  indication  that  the  transition  probabilities  are  inconsistent  and 
should  be  redefined.  This  solution  procedure  bears  some  resemblance  to 
relaxation  labeling  techniques  [10,11]  in  that  the  state  assignments  are 
iteratively  updated  by  considering  the  state  assignments  of  neighboring 
pixels.  In  practice,  the  procedure  has  been  found  to  converge  typically 
after  10  to  20  iterations. 


IV.  RESULTS  OF  SEGMENTATION 


Figure  3(a)  shows  a  digitized  aerial  photograph  of  a  rural  area  containing 
some  trees  and  fields.  The  digitized  image  is  128  by  128  pixels  in  ize 
with  gray  levels  represented  on  a  scale  of  0  to  ^55  (8  bits).  It  was 
desired  to  segment  the  image  into  two  regions  corresponding  to  the  trees 
and  fields  without  further  distinguishing  between  different  types  of  trees 
or  fields.  The  data  used  to  design  the  4  bv  4  pixel  whitening  filters  for 
each  im.ige  class  are  shown  in  the  white  boxes.  The  resulting  filters  were 
then  applied  to  the  entire  image  to  perform  the  segmentation. 

The  result  of  ML  segmentation  is  shown  in  Figure  3(b)  where  points  in 
the  image  that  were  assigned  to  tree  regions  are  c, ded  as  black  and  points 
assigned  to  field  regions  are  coded  as  white.  Although  the  true  tree  and 
field  regions  are  perceptually  descemible,  the  result  is  very  "spotty", 
in  fact,  the  Ml.  estimate  consists  of  not  two  but  a  multiplicity  of  smaller 
regions.  This  result  is  not  surprising  based  on  our  earlier  observations 
about  the  MI.  estimate. 

Figure  1(c)  shows  the  MAI’  segmentation  of  the  same  scene.  For 
this  ami  the  remaining  examples,  the  transition  probabilities  were 
taken  to  be  proportional  to  the  number  of  black  or  white  pixel  assign¬ 
ments  iii  a  square  region  surrounding  the  given  pixel.  In  this  case,  the 
state  support  region  was  taken  to  be  9  bv  9  pixels  and  the  number  of 
iterations  used  was  16.  The  MAP  result  clearly  shows  the  effect  of  the 
transition  probabilities  in  providing  connectivity  in  the  region  estimates. 
There  are  only  a  relatively  small  number  of  isolated  mi  sc  1  ass  1  tied  points. 

For  practical  purposes,  these  points  can  be  eliminated  and  the  area  near  the 
boundary  can  be  made  smooth  bv  post  filtering.  This  will  be  illustrated  in 
.i  later  example. 

Figure  4(a)  shows  another  image  consisting  of  regions  of  fields  and 
trees.  Figure  4(b)  shows  the  MAP  segmentation  of  this  image  using  the  filters 
designed  for  the  scene  of  Figure  3(a).  It  is  not  surprising  that  the  result 
Is  not  perfect  because  there  are  some  significant  differences  in  the  gray  tones 


Fig.  3.  Segmentation  of  an  image  Into  regions  of  trees  and  fields 
(a)  Digitized  aerial  photograph  (128  by  128  pixels)  showing  train¬ 
ing  data  (b)  ML  segmentation  (c)  MAP  segmentation. 
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Fig.  4.  Segment. it  ion  of  another  tree  and  field  image  vising  filter 
designed  on  data  Fig.  'l(a)  :  (n>  Digitized  photograph  (128  bv  128 
pixels)  (b)  MAP  segmentation. 


Fig.  ’  •  Steps  in  segment.it  inn  ot  two  images  of  trees  and  I  ields 
fa)  Original  digitized  images  (10  I’reproeessed  images  shew  train 
ing  data  for  MAI’  segmental  ion. 


(.1) 


Fig.  5.  Steps  in  segmentation  of  two  images  of  trees  and  fields: 
(c)  MAP  segmentation  of  preprocessed  images  (d)  MAP  segmentation 
after  post  median  filtering. 


of  each  class  between  the  two  original  images.  The  problem  can  be  alleviated 
by  preprocessing  both  images  prior  to  segmentation  to  render  the  average 
tonal  value  in  local  regions  a  middle  gray.  The  preprocessing  involves 
median  filtering  to  minimize  artifacts  in  the  boundary  regions  between 
terrain  segments  followed  by  mean  filtering  to  insure  that  the  result  will 
be  close  to  a  middle  gray.  In  particular,  the  preprocessed  image  H  is 
defined  by  the  equations 

H'(n,m)  =  F(n,m)  -  Median{F(n ,m) } 


A  (n,m)  =  ~  £  H'  (n,m) 


T  +1 

H(n,m)  =  H*  (n,m)  -A(n,m)  +  — ~2S — 


where  &  is  a  small  region  surrounding  the  point  (n,m),  M  is  the  number  of 
points  in  and  Tmax  is  the  largest  tonal  value  in  the  image  representation, 
i.e.,  that  corresponding  to  pure  white.  For  these  experiments,  was  taken 
to  be  a  square  region  11  by  11  pixels  in  size  and  Tmax  was  equal  to  255. 

This  preprocessing  is  a  modified  version  of  a  more  general  procedure  that 
can  be  used  in  various  applications  to  locally  expand  or  compress  image 
contrast  [12], 

Figure  5  shows  the  steps  involved  in  segmentation  of  the  two  images. 

Figure  5(a)  shows  the  two  original  images.  Figure  5(b)  shows  the  pre- 
processed  images  with  the  data  used  for  filter  design  shown  in  the  white 
boxes.  Figure  5(c)  shows  the  MAP  segmentation  results  and  Figure  5(d)  shows 
the  final  segmentation  results  after  post  median  filtering.  Except  for  a 
small  area  on  the  left  border  of  the  left  hand  image  (probably  due  to  some 
local  image  abnormality  not  compensated  for  by  the  preprocessing)  the  final 
results  show  a  clean  and  accurate  segmentation  of  the  images. 


Another  experiment  involved  further  segmentation  within  one  of  the 
segments  determined  in  the  previous  experiments.  The  area  of  consideration 
is  shown  in  the  white  box  in  Figure  6(a)  and  enlarged  in  Figure  6(b).  Within 
the  tree  region,  two  types  of  textures  can  be  discerned  which  correspond  to 
small  and  large  trees.  Figure  7(a)  shows  the  MAP  segmentation  of  the  tree 
area  where  the  large  trees  are  coded  in  black  and  the  remainder  of  the  image 
is  coded  in  white.  The  boundary  is  outlined  in  gray  and  superimposed  in 
white  on  the  original  image  in  Figure  7(b).  The  result  was  found  to  be  very 
accurate  when  compared  to  human  perception  of  the  tree  boundary  on  the 
original  photographic  data. 

At  the  time  of  this  report,  experiments  involving  more  than  two  classes 
of  terrain  have  not  been  performed.  However,  the  two-class  experiments 
described  suggest  that  segmentation  might  well  be  carried  out  in  a  hierarchical 
layered  form  where  an  initial  segmentation  into  gross  categories  is  first 
performed  and  then  addition  segmentation  is  performed  within  the  segments 
as  necessary.  Further  experiments  will  have  to  be  conducted  to  determine 
the  trade-offs  between  such  a  layered  approach  and  the  direct  segmentation 
into  a  number  of  classes  in  a  single  step. 
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Fig.  6.  Photograph  used  for  tree  segmentation  experiment: 
(a)  Original  digitized  image  (b)  Enlarged  section  of  image 
used  for  experiment. 


Fig.  7.  MAP  segmentation  for  tree  region  of  Fig.  6(b): 
(a)  MAP  segementation  (b)  Image  with  superimposed  tree 
boundary . 


V.  CONCLUSIONS 


The  segmentation  of  arbitrary  images  is  a  difficult  problem  because  of 
the  wide  diversity  of  image  data  and  the  corresponding  difficulty  to  estab¬ 
lish  models  of  sufficient  generality.  When  the  class  of  images  can  be 
restricted,  model  development  is  more  feasible  and  segmentation  algorithms 
can  be  designed  that  have  a  firm  theoretical  foundation.  Images  of  natural 
terrain  form  such  a  class  and  their  segmentation  is  an  important  aspect  of 
automated  image  analysis  for  reconnaissance  and  surveillance,  cartography 
and  other  similar  applications. 

This  report  described  a  class  of  models  for  terrain  images  with  two 
levels  of  structure.  The  underlying  structure  is  based  on  stochastic 
filtering  concepts  and  provides  a  means  of  representing  the  texture-like 
quality  in  local  regions  of  terrain.  Superimposed  on  this  structure  is  a 
Markov  random  field  that  describes  transitions  from  one  region  type  to  another 

Using  the  models  developed  here,  we  considered  segmentation  as  a  region 
estimation  problem  and  explored  maximum  likelihood  and  maximum  a  posteriori 
estimation  procedures.  The  ML  approach  ignores  the  Markov  structure  that 
describes  the  occurrence  of  regions  and  yields  an  algorithm  which  is  simple 
to  apply  but  produces  "spotty"  results.  The  MAP  approach  leads  to  a  set  of 
conditions  that  must  be  solved  by  iteration  but  produces  a  more  connected 
and  generally  more  desirable  result.  Examples  of  both  procedures  were  given 
on  aerial  photographic  data  of  natural  terrain.  The  examples  suggested  that 
segmentation  could  be  effectively  carried  out  with  a  layered  strategy 
beginning  with  a  gross  segmentation  of  the  image  and  proceeding  to  more 
detailed  levels  of  segmentation  within  the  original  segments.  Further 
research  would  be  necessary  to  study  the  trade-offs  involved  in  this 
approach  versus  an  approach  that  segments  the  image  into  many  different 
categories  in  a  single  step. 
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